MODELOS LINEALES

Los modelos lineales son ampliamente utilizados en acuicultura para explicar, modelar o predecir la relación lineal de una variable respuesta \(Y\) con una o mĆ”s variables predictoras \(X_1,…,X_p\).

\[Y_{i} = \beta_{0} + \beta_{1} X_{i1} + \beta_{2} \beta_{i2} + ... + \beta_{p} X_{ip} + \epsilon_{i}\]

Supuestos: \(Y_i\) se distribuye normalmente

  • Si p = 1, el modelo es una regresión lineal simple.

  • Si p > 1, el modelo es una regresión lineal mĆŗltiple.

  • Si p > 1 y alguna variable predictora es categórica, el modelo se denomina ANCOVA.

PRUEBA DE HIPƓTESIS: COEFICIENTE DE REGRESIƓN E INTERCEPTO

Prueba de hipótesis del coeficiente de regresión y el intercepto

La hipótesis nula en ambos casos es que el coeficiente del intercepto (\(β_0\)) y de la pendiente (\(β_1\)) son iguales a 0.

\(H_0:β_0 = 0\) y \(H_0:β_1 = 0\)

PREDICCIƓN A PARTIR DEL MODELO DE REGRESIƓN LINEAL SIMPLE

El Coeficiente de determinación (R²) se calcula como el cuadrado del coeficiente de correlación de pearson. Este nos indica que tan buena es la predicción que el modelo hace de los datos.

El R cuadrado ajustado (R² ajustado) nos dice qué porcentaje de la variación de la variable dependiente es explicado por la o las variables independientes de manera conjunta. En un modelo de regresión lineal, el R² ajustado es una medida de bondad que considera el número de variables existentes en el modelo.

\[R^2_{ajust} =1-(1-R^2)\frac{n-1}{n-p-1}\]

donde:

\(n\) = tamaƱo de la muestra

\(p\) = cantidad de variables predictoras en el modelo

Prueba de hipótesis del modelo completo

La hipótesis nula es si los coeficientes son iguales a 0.

\(H_0:β_j = 0\) ; \(j = 1, 2,...,k\)

Objetivos de aprendizaje

Los objetivos de aprendizaje de esta guĆ­a son:

1. - Realizar anƔlisis de datos usando modelos lineales.

2. - Realizar grƔficas avanzadas con ggplot2.

3. - Elaborar un reporte dinƔmico en formato html con Rmarkdown.

EJERCICIOS

ESTUDIO DE CASO: MUTACIƓN GEN MIOSTATINA

En vertebrados, la miostatina (MSTN) actua como es un represor del crecimiento muscular. Muchos trabajos han descrito que mutacioones en ese gen producen un fenotipo llamado de doble músculo caracterizado por un crecimiento exacerbado del músculo esquelético. En este ejercicio trabajaremos con un set de datos simulado de perros de carrera que poseen este distintivo fenotipo Adaptado de Mosher et al. 2007

Figura 1: Comparación de animales con diferentes genotipos del gen Miostatina: (A) +/+; (B) 1 mutante con alelo cys → stop (mh/+); (C) Dos copias de la mutación cys → stop (mh/mh).

Figura 1: Comparación de animales con diferentes genotipos del gen Miostatina: (A) +/+; (B) 1 mutante con alelo cys → stop (mh/+); (C) Dos copias de la mutación cys → stop (mh/mh).

La variable respuesta se denomina mass-to-height ratio y se expresa en Kg/cm. Para este set de datos existen dos variables predictoras que se deberƔn analizar: 1) Geno: que corresponde al genotipo de los animales; 2) Sex: que corresponde al sexo.

La variable Geno tiene 3 niveles denominados 0, 1 y 2 equivalentes a los genotipos +/+, +/mh y mh/mh donde 0, 1 y 2 representan el número de alelos mh. Esta forma numérica de representar los genotipos es conveniente para analizar la variable Geno como una variable cuantitativa discreta 0, 1 y 2 o como un factor (ej. factor(Geno)) según sea nuestro interés. La variable Sex tiene dos niveles male y female.

Ejercicio 1. Elaborar y configurar Reporte en formato .Rmd

Elabore un documento .Rmd y configure su reporte para exportar en .html. Instale los paquetes readxl, ggplot2, dplyr y multcomp para el anƔlisis de los datos.

knitr::opts_chunk$set(echo = TRUE)
library(readxl)
library(ggplot2)
library(dplyr)
library(multcomp)

Ejecute cada uno de los siguientes ejercicios en uno o mÔs bloques de códigos diferentes. Sea ordenado y documente su reporte adecuadamente.

Ejercicio 2. Exploratorio set de datos snp.data

Importe el set de datos myostatin.deletion.data.xlsx y transforme solo la variable Sex a factor, luego realice un anƔlisis exploratorio de datos.

Incluya:

a). Resumen estadĆ­stico de todas las variables.

snp.data <- read_excel("myostatin.deletion.data.xlsx")

snp.data$Sex <- as.factor(snp.data$Sex)
summary(snp.data)
##      animal         Geno            Sex      mass-to-height ratio
##  Min.   :  1   Min.   :0.0000   female:494   Min.   :0.2009      
##  1st Qu.:174   1st Qu.:0.0000   male  :199   1st Qu.:0.2919      
##  Median :347   Median :1.0000                Median :0.3401      
##  Mean   :347   Mean   :0.9769                Mean   :0.3590      
##  3rd Qu.:520   3rd Qu.:2.0000                3rd Qu.:0.4100      
##  Max.   :693   Max.   :2.0000                Max.   :0.5993

b). Tabla de frecuencia combinada Geno por Sex, esto es importante para ver si los datos son o no balanceados.

table(snp.data$Sex, snp.data$Geno)
##         
##            0   1   2
##   female 155  80 259
##   male   120  79   0
tapply(snp.data$`mass-to-height ratio`, factor(snp.data$Geno), mean)
##         0         1         2 
## 0.2829893 0.3263597 0.4597127
tapply(snp.data$`mass-to-height ratio`, factor(snp.data$Sex), mean)
##    female      male 
## 0.3756034 0.3177426

c). Histograma de la variable respuesta mass-to-height ratio.

ggplot(snp.data, aes(x=`mass-to-height ratio`))+
  geom_histogram(color="darkblue", fill="lightblue", bins = 12)

Ejercicio 3. Regresión lineal simple.

a). Realice una grÔfica de dispersión (geom_point) de weight en función de Geno, incluya el comando geom_smooth(method=lm) para agregar la línea de regresión a la grÔfica.

p <- ggplot(snp.data, aes(x = Geno, y = `mass-to-height ratio`))
p + geom_point() + xlab("Reference allele count") + geom_smooth(method=lm)

b). Realice un anÔlisis de regresión lineal para investigar la asociación entre weight y Geno usando las funciones lm(), summary(). Interprete los resultados del modelo lineal y responda las siguientes preguntas. ¿Qué representa el estimador de Geno y del intercepto?. ¿La pendiente es significativamente distinta de cero?. ¿Y el intercepto?

lm.geno <- lm(`mass-to-height ratio` ~ Geno, data = snp.data)
summary(lm.geno)
## 
## Call:
## lm(formula = `mass-to-height ratio` ~ Geno, data = snp.data)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.110802 -0.044236 -0.006808  0.044296  0.150224 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 0.272969   0.003434   79.49   <2e-16 ***
## Geno        0.088052   0.002615   33.67   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.06041 on 691 degrees of freedom
## Multiple R-squared:  0.6213, Adjusted R-squared:  0.6208 
## F-statistic:  1134 on 1 and 691 DF,  p-value: < 2.2e-16

c). Realice la predicción de la media poblacional de weight para los genotipos AA, AB y BB y la predicción del peso weight de tres nuevos individuos con genotipo AA, AB y BB. ¿Qué tan buena son las predicciones obtenidas? Analice el R² ajustado y los intervalos de confianza.

predict.lm(lm.geno, newdata=data.frame(Geno=c(0,1,2)), interval="confidence")
##         fit       lwr       upr
## 1 0.2729690 0.2662268 0.2797112
## 2 0.3610212 0.3565142 0.3655282
## 3 0.4490734 0.4421530 0.4559937
predict.lm(lm.geno, newdata=data.frame(Geno=c(0,1,2)), interval="prediction")
##         fit       lwr       upr
## 1 0.2729690 0.1541725 0.3917655
## 2 0.3610212 0.2423306 0.4797118
## 3 0.4490734 0.3302666 0.5678801

d). Realice una anova del modelo lineal con la función anova() (supone varianzas homogeneas) y compare con una anova que no supone varianzas homogeneas oneway.test()

# anova suponiendo varianzas iguales
anova(lm.geno)
## Analysis of Variance Table
## 
## Response: mass-to-height ratio
##            Df Sum Sq Mean Sq F value    Pr(>F)    
## Geno        1 4.1373  4.1373  1133.8 < 2.2e-16 ***
## Residuals 691 2.5215  0.0036                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# anova suponiendo varianzas diferentes
oneway.test(`mass-to-height ratio` ~ Geno, data = snp.data)
## 
##  One-way analysis of means (not assuming equal variances)
## 
## data:  `mass-to-height ratio` and Geno
## F = 520.26, num df = 2.00, denom df = 427.93, p-value < 2.2e-16

e) El siguiente modelo lineal permite realizar pruebas de contraste entre los tres diferentes genotipos.

lm(weight ~ -1 + factor(Geno), data = snp.data)

Compare los Betas y el \[R^2\] de este modelo con el realizado en el ejercicio b.

lm.geno2 <- lm(`mass-to-height ratio` ~ -1 + factor(Geno), data = snp.data)

summary(lm.geno2)
## 
## Call:
## lm(formula = `mass-to-height ratio` ~ -1 + factor(Geno), data = snp.data)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.119005 -0.044944 -0.000312  0.040813  0.139584 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## factor(Geno)0 0.282989   0.003461   81.75   <2e-16 ***
## factor(Geno)1 0.326360   0.004552   71.69   <2e-16 ***
## factor(Geno)2 0.459713   0.003567  128.89   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.0574 on 690 degrees of freedom
## Multiple R-squared:  0.9763, Adjusted R-squared:  0.9762 
## F-statistic:  9478 on 3 and 690 DF,  p-value: < 2.2e-16
summary(lm.geno)
## 
## Call:
## lm(formula = `mass-to-height ratio` ~ Geno, data = snp.data)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.110802 -0.044236 -0.006808  0.044296  0.150224 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 0.272969   0.003434   79.49   <2e-16 ***
## Geno        0.088052   0.002615   33.67   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.06041 on 691 degrees of freedom
## Multiple R-squared:  0.6213, Adjusted R-squared:  0.6208 
## F-statistic:  1134 on 1 and 691 DF,  p-value: < 2.2e-16

f) Investigue y aplique los comandos contrMat() y glht para realizar una prueba de contrastes con ajuste de bonferroni para comparaciones mĆŗltiples.

# Elabora matriz de contrastes para el factor Geno
contrastes = contrMat(table(snp.data$Geno), type="Tukey")

# Realiza comparaciones multiples
mc2 = glht(lm.geno2, linfct =contrastes)
summary(mc2, test=adjusted("bonferroni"))
## 
##   Simultaneous Tests for General Linear Hypotheses
## 
## Multiple Comparisons of Means: Tukey Contrasts
## 
## 
## Fit: lm(formula = `mass-to-height ratio` ~ -1 + factor(Geno), data = snp.data)
## 
## Linear Hypotheses:
##            Estimate Std. Error t value Pr(>|t|)    
## 1 - 0 == 0 0.043370   0.005719   7.584 3.26e-13 ***
## 2 - 0 == 0 0.176723   0.004970  35.556  < 2e-16 ***
## 2 - 1 == 0 0.133353   0.005783  23.059  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- bonferroni method)

Ejercicio 4. AnƔlisis de covarianza

a). Realice una grÔfica de dispersión (geom_point) de mass-to-height ratio en función de Geno y Sex, incluya el comando geom_smooth(method=lm) para agregar la línea de regresión a la grÔfica.

q <- ggplot(snp.data, aes(x = Geno, y = `mass-to-height ratio`, shape=Sex, color=Sex))
q + geom_point() + xlab("Reference allele count") + geom_smooth(method=lm)

b). Realice un anÔlisis de covarianza para investigar la asociación entre mass-to-height ratio, Geno y Sex usando las funciones lm(), summary().

lm.geno.sex <- lm(`mass-to-height ratio` ~ Geno + Sex, data = snp.data)
summary(lm.geno.sex)
## 
## Call:
## lm(formula = `mass-to-height ratio` ~ Geno + Sex, data = snp.data)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.107264 -0.046971 -0.006527  0.041388  0.151325 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 0.264638   0.004396  60.205  < 2e-16 ***
## Geno        0.091667   0.002864  32.006  < 2e-16 ***
## Sexmale     0.016714   0.005555   3.009  0.00272 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.06006 on 690 degrees of freedom
## Multiple R-squared:  0.6262, Adjusted R-squared:  0.6251 
## F-statistic:   578 on 2 and 690 DF,  p-value: < 2.2e-16

c). Interprete los resultados del modelo lineal.

Ejercicio 5. Anova 2 vías con interacción

a). Realice grÔficas del tamaño de los efectos plot.design() y de interacción interaction.plot de Geno y Sex sobre mass-to-height ratio.

plot.design(snp.data$`mass-to-height ratio` ~ factor(snp.data$Geno) + factor(snp.data$Sex))

interaction.plot(factor(snp.data$Geno), factor(snp.data$Sex), snp.data$`mass-to-height ratio`, mean)

b). Realice un anÔlisis de varianza con interacción para investigar la asociación entre weight, Geno y Sex usando las funciones lm(), summary().

lm.qtl <- lm(`mass-to-height ratio` ~ factor(Geno) * factor(Sex), data = snp.data)
summary(lm.qtl)
## 
## Call:
## lm(formula = `mass-to-height ratio` ~ factor(Geno) * factor(Sex), 
##     data = snp.data)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.119005 -0.041228  0.000675  0.040398  0.139584 
## 
## Coefficients: (1 not defined because of singularities)
##                               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                   0.271157   0.004495  60.321  < 2e-16 ***
## factor(Geno)1                 0.034509   0.007704   4.479 8.78e-06 ***
## factor(Geno)2                 0.188556   0.005683  33.177  < 2e-16 ***
## factor(Sex)male               0.027117   0.006805   3.985 7.47e-05 ***
## factor(Geno)1:factor(Sex)male 0.014534   0.011185   1.299    0.194    
## factor(Geno)2:factor(Sex)male       NA         NA      NA       NA    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.05597 on 688 degrees of freedom
## Multiple R-squared:  0.6764, Adjusted R-squared:  0.6745 
## F-statistic: 359.5 on 4 and 688 DF,  p-value: < 2.2e-16

c). Interprete los resultados del modelo lineal.

Ejercicio 6. Evalue supuestos del modelo

a). Evalue mormalidad del modelo con interacción del ejercicio 5 usando métodos grÔficos y prueba de shapiro.

shapiro.test (residuals (lm.qtl))
## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(lm.qtl)
## W = 0.99161, p-value = 0.0005773
plot(lm.qtl, which = 2)

b). Evalue homogeneidad de varianzas usando el comando plot.

plot(lm.qtl, which = 1)